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INTRODUCTION 


In  Refs.  [1,2,3]  the  authors  have  presented  a  "moving  singular-element" 
procedure  for  the  dynamic  analysis  of  fast  crack-propagation  in  finite  bodies. 
In  this  procedure,  a  singular-element,  within  which  a  large  number  of  analyt¬ 
ical  eigen-functions  corresponding  to  a  steadily  propagating  crack  are  used 
as  basis  functions  for  displacements,  may  move  by  an  arbitrary  amount  AT  in 
each  time-increment  At  of  the  numerical  time- integration  procedure.  The  moving 
singular  element,  within  which  the  crack-tip  always  has  a  fixed  location,  re¬ 
tains  its  shape  at  all  times,  but  the  mesh  of  "regular"  (isoparametric)  finite 
elements,  surrounding  the  moving  singular  element,  deforms  accordingly.  An 
energy-consistent  variational  statement  was  developed  in  [1,2]  as  a  basis  for 
the  above  "moving  singularity-element"  method  of  fast  fracture  analysis.  It 
has  been  demonstrated  in  [1,2]  that  the  above  procedure  leads  to  a  direct  eval¬ 
uation  of  the  dynamic  K-factor  (s) ,  in  as  much  as  they  are  unknown  parameters 
in  the  assumed  basis  functions  for  the  singular-element. 

Solutions  to  a  variety  of  problems,  were  obtained  by  using  the  above  pro¬ 
cedure  and  were  discussed  in  detail  in  [1,3].  These  problems  included,  among 
others:  (i)constant-velocity,  self-similar  propagation  of  a  finite  central 

crack  in  a  finite  panel  [analogous  to  the  well-known  Broberg's  problem]  (it) 
stress-wave  loading  of  a  stationary  central  crack  in  a  finite  panel  [analogous 
to  the  well-known  problems  of  Baker;  Sih  et  al;  and  Thau  et  al],  (it i)  constant 
velocity  propagation  of  a  central  crack  in  a  panel,  wherein,  the  propagation 
starts  at  a  finite  time  after  stress-waves  from  the  loaded  edge  reach  the  crack 
and  ( iv )  constant  velocity  propagation  of  an  edge-crack  in  a  finite  panel,  whos 
edges  parallel  to  the  crack  are  subjected  to  prescribed,  time-independent,  dis¬ 
placements  in  a  direction  normal  to  the  crack-axis  analogous  to  the  well-known 
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problems  of  Nilsson  [4].  In  Ref.  [5],  the  results  of  numerical  simulation 
of  experimentally  measured  crack-tip  versus  time  history  in  a  rectangular- 
double-cantilever-beam  (RDCB) ,  as  reported  by  Kalthoff  et  al  [6],  were  re¬ 
ported.  Also,  in  Ref.  [51,  the  authors'  results  for  the  computed  dynamic  K- 
factor  for  RDCB  were  compared  with  the  experimental  (caustics)  results  of 
[6],  and  the  independent  numerical  results  of  Kobayashi,  [7]  who  uses  a  node¬ 
release  technique  in  fast  fracture  simulation. 

In  this  paper,  the  authors  wish  first  to  clarify  and  comment  upon  sev¬ 
eral  aspects  of  the  model  formulation  which  appear  in  their  Refs.  [1,2]. 

Second,  the  model's  accuracy  and  efficiency  are  evaluated  in  terms  of  less 
sophisticated  models.  Finally,  the  practicality  of  the  special  singular 
element  for  predicting  crack  growth  for  a  given  crack  growth  criterion  is 
illustrated.  In  addressing  the  second  topic,  attention  is  focused  on:  (■£) 
the  effect  of  using  only  the  stationary-eigen  functions  (or  the  well-known 
Williams'  solution)  in  the  moving  singular-element  for  dynamic  crack  propagation, 
and  (vi)  the  use  of  isoparametric  elements  with  mid-side  nodes  shifted  so  as 
to  yield  the  appropriate  (r  )  singularity  L 8 , 9 ] .  Finally,  some  recent  results 
are  presented  which  illustrate  the  facility  of  the  propagation-eigen-function 
singular  element  for  predicting  crack  propagation  behavior  based  on  K 
versus  crack  propagation  speed  as  a  crack  growth  criterion. 

The  Variational  Principle 

Since  the  details  of  the  formulation  are  presented  in  Refs.  [1,2],  only 
those  portions  of  the  formu lat ion  necessarv  to  the  present  discussion  will 
be  include  1  here.  Further,  for  simplicity  the  rather  general  equations  of 
Ret  ..  !  1 ,  1  will  be  spec  i  ilir.ed  to  the  case  of  Mode  [  crack  growth  in  bodies 
sue  1  •  s’ t  to  traction  tree  crack  .urfaccs,  .tern  body  torcc  and  with  geometry/ 


applied  loading  such  that  the  model  can  make  use  of  symmetry  about  the  crack 
plane.  In  Ref.  [2],  the  principal  of  virtu  l  work  proposed  as  the  governing 
equation  for  crack  propagation  during  the  period  which  the  crack 

elongates  by  amount  A£  is  given  as: 


0  -j  j(o  ?  ,+o  ^  ,)6e^ .  4-  p(ii2+iiW}  dV 
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The  superscripts  1  and  2  refer  to  quantities  at  times  t^  and  t  respectively. 

The  integrals  in  order  of  appearance  refer  to  the  volume  of  the  body  at  time 

t^.  the  portion  of  the  surface  of  the  body  at  time  t^  subjected  to  prescribed 

tractions  and  the  new  crack  surface  created  between  times  t^  and  t^.  Note 

2  2 

that  the  variational  quantities  Se^  and  Su^  reflect  the  kinematic  constraint 
at  t2  and  therefore  are  arbitrary  on  A£.  Making  use  of  the  small  strain  dis¬ 
placement  relation,  the  symmetry  of  o^and  using  the  divergence  theorem,  the 
first  term  of  the  volume  integral  in  (1)  becomes: 
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(2) 

where  3V2  is  the  boundary  of  Noting  that  (i)  A£  is  part  of  3V^  (resulting 

in  the  last  term  of  (1)  dropping  out),  (H)  that  is  a  part  of  3V^  and  (Hi) 

2  2 
that  5u^  is  zero  on  any  portion  of  3V  that  has  prescribed  displacements,  we  have 
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It  is  seen  that  (4c)  is  nothing  other  than  the  condition  that  the  new 

crack  surface  be  traction  free.  Equation  (4b)  stipulates  that  traction 

boundary  conditions  be  satisfied  on  and  (4a)  is  a  statement  of  dynamic 

^  1  ..1 

equilibrium  within  the  body.  If  one  assumes  that  a..  ,=pu,  then  (4a)  reduces 

300  1 

2.2  2  ..2 

to  the  usual  equilibrium  equation  at  t„  (ie,  a..  .=pii.).  Then  since  a..  ,=pu. 

2  1JO  i  1J,J  i 

it  follows  that  the  state  at  t^  must  also  satisfy  the  usual  equilibrium  equa¬ 
tion  and  so  forth.  Therefore,  (4a)  is  equivalent  to  the  standard  expression 
for  dynamic  equilibrium  when  the  assumption  that  o^  is  val3-d*  A  similar 

argument  leads  to  (4b)  reducing  to  the  usual  condition  for  satisfaction  of  trac¬ 


tion  boundary  conditions. 

In  the  finite  element  model  formulation,  the  above  assumption  is  not 

valid  and  therefore,  the  equations  (4)  do  not  reduce  to  those  usually  found  in 

finite  element  model  derivations.  One  reason  for  o^du!  is  that  in  modeling 

3.J  t 

crack  growth,  it  becomes  necessary  in  the  procedure  of  Refs.  [1,2]  to  change  the  mesh 

configuration  at  each  crack  growth  time  step  and  to  interpolate  displacement, 

velocity  and  acceleration  data  at  the  new  node  locations.  This  interpolation 

will  in  general  result  in  some  disequilibrium  in  the  interpolated  solution. 

If.  :  or  example,  we  assume  some  disequilibrium  at  t,  sucli  that  .  =f , 

1  3-  3.J.J 

■  2  2 

then  ;.iL  is  I  act  ion  of  (4a)  leads  to  ,iu‘"-i7.  ,--t.  Clearly,  the  disequilibrium 

i  i  J  ,  1 

at  aihsequent  steps  will  be  the  same  form  (f)  with  the  sign  alternating  at 
eacn  ste:>.  Thereiore,  it  would  appear  that  the  formulation  of  Ref.  [1,2]  should 
reswi  t.  in  o  :c  i  I  lat  ions  . 
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Numerical  experimentation  with  the  formulation  has  indeed  shown  this 
oscillation  to  occur.  However,  when  used  with  the  special  singular  element 
of  Refs.  [1,2]  the  only  time  it  has  occured  at  discernable  levels  is  in 
static  analyses.  In  dynamic  analyses,  it  has  been  found  that  the  inertial 
forces  are  generally  large  enough  to  make  the  oscillatory  forces  negligible. 

Similar  oscillation  has  been  observed  when  implementing  the  proposed 
principle  of  virtual  work  with  eight-noded  isoparametric  elements.  (Wherein 
the  appropriate  crack  tip  singularity  was  obtained  by  shifting  midside  nodes 
as  suggested  by  Ref.  [8,9]).  It  is  generally  found  that  the  oscillations 
when  using  the  isoparametric  elements  are  larger  than  those  observed  with  the 
special  singular  element.  This  is  believed  to  be  the  result  of  inherently 
larger  interpolation  induced  disequilibrium  with  these  less  sophisticated 
elements. 

A  related  variational  principle  for  quasi-static  crack  growth  in  elastic- 
plastic  bodies  was  also  presented  in  [1],  To  account  for  effects  of  history 
dependent  plasticity  and  finite  deformation  gradients,  an  updated  Lagrangean 
rate  formulation  was  used.  As  a  result  of  the  formulation  yielding  an  in¬ 
cremental  analysis  (as  opposed  to  the  computation  of  total  state  quantitites 
as  in  the  elasto-dynamic  analysis  above) ,  the  effect  of  state  quantities  at 
t^  appearing  in  the  finite  element  equations  for  the  solution  at  t^  is  funda¬ 
mentally  different.  The  variational  statement,  simplified  for  the  case  of 
zero  body  force,  traction  free  crack  surfaces  and  symmetrical  modeling  of  Mode 
I  crack  growth.  Is  given  as: 


f  [*r  +  S  5  IdV  -  f  T6u.dS  +f  r  1  .v  L6  u  .  d  E  =  0  (5) 

L  lJ  ij  Li  Js  1  1  1  *-J  J  l 

wnero  1 1  arc  the  Cauchy  stress  components  of  the  solution  at  t^  in  the 
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current  reference;  eoitf  igurnt  ion  (t  iV  tiro  the  incremental  displace¬ 
ments  in  going  from  the  rctereiice  state  to  the  final  state  ( t2  ’  ^2  ’  > 

T.  are  the  applied  incremental  tractions,  g..=ti  .  u  g..=u.  ,+u.  .  and 

L  1J  K.,1  K.,j  ij  i,j  j,i 

S_  are  the  incremental  second  Piola-Kirchhof f  stress  components  (or  what  are  also 
known  as  Truesdell  stress  increments).  Equation  (5),  through  the  use  of  the  di¬ 
vergence  theorem,  leads  to  the  following  Euler-Lagrange  equations: 

(rf.u.  .),.+§..  .  +  F.  =  0  in  V,  (6a) 

kj  i,k  j  ij ,j  l  1 

1  •  - 

(t,  .u.  +S..)v.  -  T.  =  0  on  S  (6b) 

kj  t,k  ij  j  l  ax 

/T-*-.u.  ,+S.  .+!?■.  )v?"  =  0  on  AE  (6c) 

(  k]  ij  ij  j 

Equation  (6a)  is  the  usual  translational  equilibrium  condition  associated 
with  updatcd-Lagrangean  formulations  in  terms  of  the  second  Piola-Kirchhof f 
stress  and  does  not  contain  any  additional  terms  such  as  found  in  equation 
(4a).  Equation  (6b)  is  the  usual  condition  for  satisfaction  of  traction 
boundary  conditions  and  equation  (6c)  is  the  condition  that  the  newly  created 
crack  surface  be  traction  free.  Further  study  of  the  derivation  of  equations 
(5)  and  (6)  shows  that  there  is  nothing  inherent  in  the  formulation  to  account 
for  (or  correct)  the  error  from  interpolation  of  quantities  from  the  finite 
element  mesh  at  t^  with  crack  length  to  the  finite  element  mesh  at  t^  with 
crack  length  E^+AE.  However,  the  accumulated  error  at  increment  P  (£p)  can 
be  measured  by: 


”f„  TP.dgL.dV  -  f  ?P  HI.  . 

Jvp  11  lJ  J  s  1  1 


IS 


wiiicu  in  a  check  on  the  equilibrium  at  increment  p.  Since  this  rate  f ormulj t ion 
lie.'  .  not  result  in  terns  which  lead  to  oscillation  of  the  solution  it  seems 
time  ier  Milat ion  of  the  1  i  near-e  las t  ic  dynamic  problem  in  terms  of  an  in- 
cre  -.ent.il  node  1  similar  to  that  for  the  e  last  ic-p  last  ic  large  deformation  problem 
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would  eliminate  the  trouble  with  oscillations  while  at  the  same  time  retain 
the  crack  growth  modeling  features. 

Results  of  analyses  which  are  presented  later  in  this  paper  are  largely 
based  on  the  formulation  as  originally  proposed  in  Ref.  [1,2]  since  it  ap¬ 
pears  that  no  significant  change  would  result  from  a  reformulation.  The  one 
exception  to  this  are  the  results  obtained  through  the  use  of  isoparametric 
elements  with  midside  nodes  shifted  so  as  give  a  singularity  at  the  crack  tip 
In  those  analyses,  the  usual  statement  of  virtual  work  is  applied: 


0  -f  ja2.6e2  +pii2SuJ  dV  -  f  T^duJ  dS 

J  I  i]  ij  i  i|  Js  i  i 


Traction  free  crack  surfaces  are  approximated  by  letting  nodal  forces  on  the 
crack  surfaces  be  zero. 

Singular  Element  for  Dynamic  Crack  Propagation 
A  singular  crack  tip  element  was  also  developed  in  [1,2]  and  used  in  con¬ 
junction  with  the  formulation  (2)  for  the  analysis  of  Mode  I  dynamic  crack 
propagation  in  linear  elastic  two  dimensional  bodies.  This  singular  element 
uses  an  arbitrary  number  of  the  displacement  eigen-functions  which  come  from 
the  solution  of  a  crack  in  an  infinite  body.  For  dynamic  crack  propagation, 
these  eigen-functions  are  taken  as  those  for  the  corresponding  steady-state 
dynamically  propagating  crack  in  an  infinite  body.  Equations  (7)  give  the 
form  for  the  assumed  displacement,  velocity  and  acceleration  within  the  sin¬ 


gular  element: 


uS(',,x2>t)  =  U(£,x2,v)S(t) 


uJ  =  L’S  -  v(U),r3 

uS  =  Vi}  -  2  v(U)  3  +  V2(U)  3 


(7a) 

(7b) 

(7c) 


where  V  is  the  matrix  of  eigen-functions  (pLus  appropriate  rigid  body 


modes ) 


is  the  vector  of  undetermined  coefficients,  v  is  the  crack  propa¬ 


gation  speed  and  are  the  coordinates  relative  to  the  moving  crack  tip 

(r=x^-vt).  It  should  be  noted  that  (7b)  and  (7c)  are  obtained  form  (7a) 
through  differentiation  with  respect  to  time  with  the  assumption  that  v  is  not 
a  function  of  time.  The  following  comments  should  remove  any  incorrect  notions 
that  this  in  any  way  limits  the  use  of  the  element  to  constant  speed  crack 
propagation.  First,  it  has  been  shown  [12]  that  the  near-tip  fields  are 
the  same  for  steady-state  and  transient  crack  propagation.  Therefore,  provided 
v  at  each  time  step  reflects  the  current  speed,  there  is  no  question  that  the 
eigen-functions  for  the  element  are  correct  and  that  the  coefficient  of  the 
singular  eigen-function  (3^)  is  indeed  the  Mode  I  stress  intensity  factor.  A 
second  cons iderat ion  is  that  the  associated  stress  eigen-functions  do  not 

satisfy  the  stress  equilibrium-equation  for  non-steady-state  (as  viewed  by  an 

observer  moving  with  the  crack-tip).  . 

32  u .  „  32u A 

3t  -  2v3^  +  v  (8) 

but  instead,  satisfy  the  corresponding  steady  equation: 


o .  .  .  =  iV 

ij  .  J 


2 

3C2 


(9) 


While  it  would  be  preferred  that  equation  (8)  be  satisfied  exactly,  the  dis¬ 
placement  finite  element  method  does  not  require  this.  The  stress  equilibrium 
of  (3)  is  satisfied  in  the  usual  approximate  sense  associated  with  the  finite 
element  method. 
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metric  elements.  This  method  involves  selecting  8  such  that  the  following 
three  error  functionals  are  minimized: 


(uS-uR)2dp;  12  =  J  (uS-uR)2dp;  I3  =  j  (oS-ii 


"S  "R)2dP  (10) 


s  s  **s  R  R  R 

where  u  ,u  ,u  are  as  in  equations  (7),  u  ,u  ii  are  displacements,  velocities, 

and  accelerations  at  the  boundary  of  the  eight-noded  isoparametric  elements, 

and  pg  is  the  boundary  of  the  singular  element  shared  by  isoparametric  elements. 

This  procedure  gives  g,  g  and  g  in  terms  of  qg,  qg  and  qg,  the  nodal  displacements, 

velocities  and  accelerations  of  nodes  on  the  boundary  shared  by  the  singular 

element  and  the  isoparametric  elements.  In  particular.  Bis  related  to  qg  by 


§  ”  b  qs 
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where  A  is  in  general,  a  rectangular  MxN  matrix  (M  being  the  dimension  of  6 
and  N  the  dimension  of  q  ) .  The  exact  form  of  A  does  not  enter  the  present 
discussion  but  can  be  found  in  [1,2].  The  question  to  be  addressed  here 
concerns  the  constraints  which  must  be  placed  on  the  dimension  of  A  (and  there¬ 
fore  3)  .  It  is  well  known  that  in  hybrid  finite  elements  there  is  a  restriction 
on  the  number  of  internal  parameters  (g)  so  that  the  matrix  equations  relating 
these  parameters  to  the  unknowns  of  the  final  finite  element  equations  are  non¬ 
singular  [11].  It  should  be  noted  that  the  current  procedure  is  not  (strictly 
speaking)  a  hybrid  procedure  and  tiierefore  the  matrix  A  relating  8  to  qg  does 
not  pose  the  same  problem.  However,  the  following  does  show  there  is  a  res¬ 
triction  on  the  dimension  of  8,  in  the  derivation  of  [1,2]  there  are  nine  nodal 
points  and  therefore  eighteen  degrees  of  freedom  associated  with  the  singular 

element.  All  of  these  are  on  its  boundary  (denoted  a  in  (10))  as  illustrated 

s 

in  the  typical  mesh  of  Figure  L.  To  ascertain  the  behavior  of  the  singular 


element  with  differing  numbers  of  e i gen- fune t ions  (ie,  differing  dimensions 


>f  the  singular  element  for  zero 


for  )  tlu>  e  igen-va  lues  and  e  i  gen-modes  i 
crack  sneed  were  calculated.  In  making  Llie  calculations,  the  symmetry  plane 
nodal  displacement  normal  to  the  crack  plane  was  constrained  since  no  eigen¬ 
functions  were  included  for  rigid  body  translation  normal  to  the  crack  plane. 
Based  on  the  deletion  of  this  one  degree  of  freedom  it  can  be  seen  that  the 
singular  element  must  have  seventeen  deformation  modes  (eigen  vectors)  and 
should  have  only  one  zero  energy  mode  (eigen  value  of  zero),  corresponding 
to  a  rigid  body  translation.  In  varying  the  dimension  of  3  (call  it  M)  it 
was.  found  that  for  M  =  17,  there  was  one  zero  eigenvalue  which  through  ex¬ 
amination  of  the  eigenvector  did  indeed  correspond  to  the.  desired  rigid  body 
translation  parallel  to  the  crack  plane.  When  cases  for  M  < 17  were  considered, 
the  number  of  these  zero  energy  mod.es.  (P)  was  P  =  18-M,  (M<_17)  .  These  excess 
zero  energy  modes  are  clearly  undesirable.  It  was  also  observed  that  the 
desired  rigid  body  mode  was  no  longer  present  when  these  extra  inodes  occured. 
The  constraint  on  >1  for  the  current  conf iguration  is  therefore  M  >  17.  In 
the  general  2D  case  where  N  Is  the  number  of  unconstrained  nodal  degrees  of 
freedom  associated  with  the  singular  element  we  must  have  M  >_  N,  where  M  in¬ 
cludes  the  appropriate  number  of  special  rigid  body  displacement  functions. 

If  we  denote  the  number  of  these  rigid  body  functions  by  R,  then  the  required 
number  of  crack  solution  eigen-functions  (C)  is  given  by  C  2.  which  is 
exact lv  that  constraint  four.'  for  the  number  of  assumed  element  stress  func¬ 
tion;  La  hvbrid  stress  finite  element  models  [111.  The  results  presented  in 
[1,  ! ,  I  1  and  to  he  presented  here  have  used  M-.’.O  therefore  satisfying  the  above 
con  . t  r a i at . 

Procedures  foty  -.jnul.it  ion  of  ’’rack  j'.rowth 
Ike  model  in;  of  .-.rack  oropag.it  ion  with  either  the  special  singular  element 
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or  with  the  singular  isoparametric  elements  requires  the  finite  element 
mesh  in  the  region  of  the  crack  tip  to  be  modified  at  each  time  step  for 
which  crack  growth  occurs.  The  procedure  used  here  (and  in  [1,3,5]  is  to 
shift  the  singular  element  (s)  without  distortion,  (as  shown  in  Figure  2)  so 
that  it  is  only  the  regular  elements  surrounding  the  special  modeling  region 
which  become  distorted.  -An  important  point  which  differentiates  this  pro¬ 
cedure  from  the  more  common  node  release  techniques  for  crack  growth  is  that 
the  size  of  the  increment  of  crack  growth  (AE)  is  not  restricted  by  nodal 
spacing  but  rather  can  be  made  as  small  as  desired.  Figure  2  also  illustrates 
that  the  distortion  of  elements  periodically  reaches  a  critical  degree  at  which 
time  the  crack  tip  region  of  the  model  is  reraeshed  before  applying  the  shifting 
procedure.  Since  this  procedure  involves  the  shifting  of  nodal  points  and 
since  the  nodal  quantities  of  displacement,  velocity  and  acceleration  at  each 
time  step  appear  in  the  difference  equations  associated  with  the  Newmark  time 
integration  scheme  for  the  solution  at  the  subsequent  time  step,  it  is  neces¬ 
sary  to  use  interpolation  procedures  to  obtain  correct  values  for  the  shifted 

4 

nodes . 

Considerations  of  Efficiency  and  Accuracy 

In  the  previous  sections,  the  discussion  was  largely  in  terms  of  partic¬ 
ular  features  of  the  variational  principle  or  singular  element  derivation. 

Here  the  discussion  will  be  of  a  broader  nature  with  most  of  the  attention 
being  focused  on  the  more  general  attributes  of  efficiency  and  accuracy. 

The  singular  element  with  crack  propagation  speed  dependent  eigen-func¬ 
tions  has  several  attractive  features  which  stem  directly  from  the  use  of  the 
analytical  solution  for  the  near-c rack-t ip  field.  First,  the  traction  free 
crack  surface  conditions  are  satisfied  exactly.  Second, the  coefficient  of  the 
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singular  eigen- fund  ion  is  the  Mode  1  stress  intensity  factor  and  is  obtained 
directly  without  recourse  to  indirect  energy  based  procedures  or  extrapolation 


methods.  Third,  because  many  eigen-functions  are  used,  the  accuracy  of  the 
.solution  is  less  sensitive  to  the  singular  element  size  than  for  elements  with 
less  elegant  basis  functions. 

As  might  be  expected,  this  element  also  has  some  features  which  tend  to 
offset  the  above  positive  ones.  One  such  feature  is  that  the  propagation- 
eigen-functions  lead  to  a  non-symmetric  stiffness  matrix.  However,  since  the 
non-symmetry  is  localized  to  rows  and  columns  corresponding  to  a^,  the  ad¬ 
ditional  effort  in  solving  the  equations  is  not  excessive.  Another  aspect  of 
using  special  elements  which  requires  attention  is  the  inherent  incompatibility 
of  displacement  with  neighboring  elements.  The  procedure  proposed  in  [1,2] 
and  used  here  involves  satisfying  compatibility  along  p^  in  an  integrated  least 
square  sense  as  shown  in  equation  (10 )  . 

In  order  to  weigh  the  above  positive  features  against  the  negative,  two 
alternative  models  are  considered.  The  first  alternative  is  to  use  a  similar 
special  singular  element  but  to  substitute  stationary  crack  eigen-functions 
for  the  propagating  crack  eigen-functions.  This  substitution  has  two  effects. 
First,  the  stiffness  matrix  becomes  symmetric.  Second,  the  coefficient  of  the 
singular  eigen-function  can  no  longer  be  interpreted  as  K^.,  the  Mode  I  stress 
intensity  factor.  At  first,  this  loss  seems  a  dear  price  to  pay  since  it  would 
appear  that  one  must  resort  to  indirect  procedures  for  determining  such 
as  energy  calculations  or  the  fitting  of  propagation-eigen-functions  to  the 
near  tielJ  solution!  However,  can  be  obtained  direetlv  from  the  s tat ionarv- 
e  i  gen- func  t  i  on  element  solution  us  in-.;  a  simple  formula  described  later.  Witu 
tni  question  of  calculating  answered,  it  remains  to  be  seen  how  well  the 
chosen  number  ot  scat  ion. irv  functions  (twentv)  can  accomodate  thi.  Jistored 
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near  field  displacement  patterns  of  a  propagating  crack. 


It  is  clear  that  whether  one  uses  stationary-eigen-functions  or  propagation- 
eigen-functions,  the  special  element  will  still  require  additional  work  to 
ensure  compatibility  with  neighboring  elements.  To  ascertain  the  benefits  of 
eliminating  this  additional  work,  a  third  model  is  considered.  This  model,  as 
mentioned  previously,  uses  eight-noded  isoparametric  elements  exclusively. 

-h 

The  r  singularity  in  stress  and  strain  is  incorporated  in  the  model  by  shifting 
midside  nodes  on  element  edges  joining  the  crack  tip  node  to  the  quarter-point 
of  the  element  side  as  illustrated  in  Figure  1.  To  ensure  the  correct  behavior 
in  all  angular  directions  from  the  crack  tip,  the  elements  adjoining  the  crack 
were  degenerated  to  the  triangular  form  seen  in  Figure  1.  While  this  model 
has  no  problem  in  terms  of  compatibility,  it  does  have  a  problem  in  terms  of 
calculation  of  K^.  One  has  virtually  no  choice  but  to  resort  to  energy  methods, 
fitting  of  eigen-functions  or  extrapolation  procedures.  Since  a  major  consider¬ 
ation  in  comparing  the  above  models  is  the  accuracy  and  ease  of  computing  K^, 
the  indirect  procedures  for  computing  with  the  stationary  eigen-function 
model  and  quarter-point  isoparametric  element  model  will  be  described. 

Indirect  Methods  for  Computing  K^. 

As  pointed  out  previously,  when  the  propagation-eigen-function  singular 
element  is  used,  is  evaluated  directly  during  the  solution  procedure  and 
therefore  indirect  methods  for  evaluation  K^.  are  not  required.  However,  that 
does  not  mean  that  these  methods  can  not  be  used  and  in  fact  the  comparison 
of  from  alternate  procedures  is  a  good  manner  for  checking  the  consistency 
of  the  solution. 

The  first  two  indirect  procedures  rely  on  the  well  known  relationship 

* 

between  and  the  fracture  energy  release  rate  ((',)  for  pure  Mode  I  fracture: 


where  f(  a  ,oe  )  =  a  (l-a2)/[4a  u  -  (l+u“")2] 
as  as  as  s 

a2  =  1  -  (v/C  )2;  a2  =  1  -  (v/C  )2 
d  a  s  s 

C2  =  u/p 
s 

C  ,  =  —  — — )  for  plane  strain 

d  p  l-dv 

or  C  ,  =  —  (g - )  for  plane  stress 

d  p  1-v 

In  the  above,  v  is  the  crack  propagation  speed,  p  the  shear  modulus,  v  the 

possons  ratio,  p  the  mass  density,  C  the  shear  wave  speed  and  C,  the  dila- 

s  d 

tational  wave  speed.  For  the  limiting  case  as  v  goes  to  zero,  we  have 

f(l,l)  =  1-v  for  plane  strain 


or  f(l,l) 


for  plane  stress 


To  make  use  of  (12),  one  must  evaluate  0.  This  has  been  done  in  the  present 
work  using  three  different  approaches.  The  first  approach  to  determining 
G  is  through  an  energy  balance.  This  is  done  in  crack  propagation  analysis 
most  easily  by  considering  the  energy  of  the  system  at  two  adjacent  time  steps 
t,  and  t  ,  between  which  an  increment  in  crack  length  A2  has  occured  (A2=2 . 
It  wo  define  the  increment  in  work  done  on  the  system  by  the  applied  traction 
and.  d  i  ..p  lavement  boundary  conditions  as  AW  =  W^-W  and  similarly  denote  the 
change  in  strain  energy  by  AU  and  the  change  in  kinetic  energy  by  AT  then  an 
average  G  for  the  interval  (t  ,t7)  is  given  hv: 

_  AW- 'T-  ■  t; 

In  ■:.}  (11) 

where  b  is  the  Length  ot  tin'  propagnt  ing.  c rack  front. 

Th.-  second  approach  to  computin'.’.  G  is  ro  use  a  crack  closure  integral. 


Vii  i  ;  calculation  involves  the  tractions  on  M  existing  at  t^  .  For  pure  Mode 


Li 


fracture  behavior  one  has  only  displacement  component  uy  and  traction  com¬ 
ponent  Ty  present  at  AE.  An  average  G  for  the  interval  t^,  t ^  is  given  by: 

G  =  ^  f  Ty(x,t1)uy(x,t2)d:c  (15) 

AE 


where  u  is  half  the  total  crack  opening  due  to  the  use  of  symmetry  in  the 
model.  The  factor  of  usually  observed  in  linear  elastic  work  evaluations 
is  canceled  by  a  factor  of  2  which  accounts  for  the  use  of  symmetry  in  the 
model . 


The  third  method  for  evaluating  G  is  the  J-integral.  It  is  well  known 
that  G=J  for  elastic  bodies  and  therefore  is  obtained  from  the  definition 
of  J  in  the  current  symmetrical  dynamic  analysis  of  a  Mode  I  crack  [10]: 


J=2 


(16) 


where  p  is  the  mass  density,  W  is  the  strain  energy  density,  T  is  a  curve 
connecting  the  upper  crack  surface  to  the  symmetry  plane  ahead  of  the  crack 
(which  is  propagating  in  the  x  coordinate  direction),  A  is  the  two-dimensional 
region  enclosed  by  i,  n^  is  the  x-eomponent  of  the  outward  unit  normal  to  T 
and  T  is  tile  traction  vector  acting  on  T (outward  positive).  The  factor  of 
2  again  reflects  'he  use  of  symmetry  in  modeling. 

In  addition  to  the  above  energy  related  procedures,  two  additional  methods 
are  used  here  for  determining  K^.  The  first  of  tlie.se  two  involves  the  fitting 
of  near  field  displacements  with  the  propagation-eigen-functions  discussed 
previously.  This  method  is  quite  general  in  nature  and  can  be  used  with  either 
thi'  special  stat i onary-e i gen-f unc t ion  element  or  with  the  quarter-point  singular 
element  ;.  The  procedure  is  to  use  equ.it  ions  (10)  to  obtain  an  equation  of  the 


for . i  i  I  1 )  such  that  the  coefficient  ot  the  singular  propagat ion-eigen-function 


3  (3  =K.  )  can  bo  related  to  the  near  field  nodal  displacements  q  . 
Ill  -s 


The  final  procedure  for  obtaining  K  is  primarily  applicable  to  the 

stationary-eigen-function  element  and  results  in  the  evaluation  of  from 

this  element  being  nearly  as  direct  as  when  using  the  propagation-eigen- 

p 

function  elemenc.  Let  6^  be  the  undetermined  coefficient  of  the  singular 

P  S 

propagation-eigen-function  (3^=K^.)  •  Let  3^  be  the  corresponding  coefficient 

S  P 

of  the  s tationary-eigen-f unctions ,  noting  3^  =  6^  =  K^.  only  if  v,  the  crack 
propagation  speed,  is  zero.  Using  the  definition  of  G  for  Mode  1  fracture: 


AT 


G  -  limit 
AZ-HD 


A£  / 


T  (x)  u  (x-AZ)  dx 

y  y 


(17) 


it  is  possible  to  obtain  two  equivalent  definitions  of  G  in  terms  of  3^  or 

S 

3^  by  substitution  of  the  respective  singular  eigen-functions  into  (17): 


c  5  liiid  (sy  and  c  5  („*)* 


(18) 


From  (18)  it  is  then  seen  that  6^  which  is  always  equal  to  can  be 


related  to  3^  by: 


K  =  3^  = 

1  1  V  ("d’as) 


(19) 


where  f(L,l)  and  fOa^.a^)  are.  as  defined  in  the  text  immediately  following 
equat ion (12 ) • 

Crack  Propagation  Comput.it  ions  for  Comparison  of  Models 
The  primary  purpose  of  these  computations  is  to  compare  results  using 
several  levels  of  modeling  sophistication.  To  keep  the  problem  from  ob¬ 
scurin',  the  basic  differences  in  modi' Ling  techniques,  a  rather  simple  com- 
hinai'.on  of  geometrv  and  loading  was  selected.  The  problem  which  was 
cons  id.  previously  in  Met'.  (  1 1  i  .  that  of  t  lie  constant  velocity  prop- 


'at  ;  on  •>:  an  edge  crack  in  a  square  sheet  whose  edges 
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parallel  to  Che  direction  of  crack  propagation  are  subjected  to  uniform 
displacements  in  the  direction  normal  to  that  of  crack  propagation. 

The  dimensions  and  mesh  configuration  are  depicted  in  Figure?  1.  This  prob¬ 
lem  is  analogous  to  that  treated  by  Nilsson  [4]  who  obtained  an  analytical 
solution  for  steady-state  stress-intensity  factor  for  the  constant  velocity 
propagation  of  a  semi-infinite  crack  in  a  finite-height,  infinite-width 

4  2 

strip.  In  the  following,  the  material  properties  are  v=0.286,  y=2.94  x  10  N/ram 
—6  3 

and  p=2.45  x  10  kg/mm  .  Three  levels  of  constant  crack  velocity  propagation, 

(v/C  )=0.2,0.4,0.6  are  considered.  In  each  case,  the  initial  crack  length, 
s 

(E  Al)  is  0.2.  For  each  case,  the  crack  growth  increment  size  is  maintained 
o 

constant  (as  opposed  to  time  step  size)  at  a  value  of  AE/W=0.005  which  cor¬ 
responds  to  2.5  percent  of  the  singular  element’s  dimension  in  the  crack  prop¬ 
agation  direction. 

In  these  analyses  the  uniform  prescribed  displacement  is  applied  sta¬ 
tically  at  t=0  sec.  While  maintaining  this  prescribed  displacement,  the  crack 
is  then  made  to  propagate  at  a  uniform  speed.  This  procedure  results  in  crack 
acceleration  over  the  first  time  step  and  since  the  time  step  size  is  varied 
so  as  to  maintain  the  crack  growth  increment  size  constant,  this  acceleration 
will  also  differ  for  each  value  of  the  uniform  crack  speed. 

v=0. 2C.  The  computated  K^(t)  for  this  lowest  considered  propagation 
speed  are  illustrated  in  Figure  3.  Note  that  the  computed  K^(t)  are  normalized 

CO 

with  respect  to  K  ,  which  is  the  static  K  for  the  infinite  widtli  strip  problem  of 

GO  _  —  2  c 

Nilsson  (whore  K^=(  u.^E)  /  *  li  ( 1-  v  )  and  are  plotted  as  a  function  nr  nondiniens  ionnl 
crack  Length  (!!/', v').  In  each  of  Figures  3,4  and  5  the  arrows  indicate  the 

crack  lengths  at  which  the  remesh i ng  procedure  (illustrated  in  Figure  2)  takes 
place.  The  dashed  Lines  of  Figures  I,  4  and  5  indicate  the  steadv-state  value 
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of  for  the  infinite  strip  problem  due  to  Nilsson  [4],  The  nearness  to 
unitv  of  the  dashed  line  of  Figure  3  indicates  the  relatively  low  velocity 
dependence  of  the  strip  solution  at  this  crack  speed. 

The  solid  curve  of  Figure  3  represents  the  values  of  fron  the  prop- 
agat ion-eigen-function  element  and  as  discussed  previously  are  determined 
directly  as  the  coefficient  of  the  singular  eigen-function  during  the  solution. 
The  solid  points  of  Figure  3  are  the  results  from  the  stationary-eigen-function 
element.  To  illustrate  the  small  effect  of  crack  speed  at  v=0.2Cs,  the  plotted 
values  are  just  the  coefficient  of  the  singular  eigen-function  but  because  of 
the  non-zero  crack  speed  are  not  strictly  speaking  values  of  .  To  obtain 
correct  values  of  K^.  in  this  case,  one  would  have  to  use  equation  (19).  Since 
this  correction  would  uniformly  lower  all  these  points  by  only  1.4%,  these 
corrected  values  are  not  plotted.  It  can  be  seen,  however,  that  this  correction 
would  indeed  tend  to  improve  the  already  excellent  agreement  with  the  propa- 
gat  ion-e i gen- function  element. 

The  open  sympols  of  Figure  3  are  the  results  of  the  quarter-point  isopara¬ 
metric  element  computations  The  open  circular  points  are  based  on  a  G  obtained 
through  a  global  energy  balance  (14)  which  is  then  converted  to  K^.  through 
equation  (13).  It  can  he  seen  that  except  for  the  first  four  time  steps,  these 
point:;  agree  quite  well  with  the  propagat ion-eigeu-funct ion  element  solution. 

The  open  triangular  points  are  based  on  a  C  obtained  through  the  crack  closure 
integral  of  equation  (Lb)  and  then  converted  to  through  equation  (12).  It  is 


clear  from  Figure  3  that  the  crack  closure  integral  procedure  is  inferior 
to  the  energy  balance  procedure  when  used  with  the  isoparametric  elements. 

This  is  attributed  to  the  inherent  inaccuracy  of  the  boundary  tractions 
for  the  isoparametric  element  which  are  required  when  using  equation  (15) • It 
should  be  noted  that  a  very  satisfactory  aspect  of  the  crack  growth  modeling 
procedure  used  here  is  that  none  of  the  curves  in  Figure  3  show  erratic  be¬ 
havior  which  can  be  related  to  the  remeshing  procedure.  From  Figure  3,  it 
seems  that  any  one  of  the  three  models  is  adequate  for  the  problem  when  v=0.2Cs- 
y=0. The  results  for  this  intermediate  case  are  given  in  Figure  4. 
Again,  the  dashed  curve  represents  the  steady-state  infinite  strip  solution. 

To  illustrate  the  increasing  effect  of  the  crack  speed  on  the  crack  tip  field, 
the  solid  points  are  again  the  coefficient  of  the  singular  term  of  the  sta- 
tionary-eigen-function  solution  and  must  be  modified  via  equation  19  before 
being  interpretable  as  .  The  solid  square  points  are  these  modified  values. 
It  can  be  seen  that  the  agreement  between  the  propagation-eigen-function 
solution  and  the  stationary  eigen-function  solution  is  still  quite  good  despite 
the  increased  effect  of  crack  speed  on  the  crack  tip  field. 

The  quarter-point  isoparametric  element  results  are  indicated  by  open 
circular  points.  Here  again,  the  global  energy  balance  procedure  has  been 
used.  Unlike  the  results  for  v=0.2C  ,  there  is  a  pronounced  difference  be¬ 
tween  the  isoparametric  element  results  and  the  special  element  results.  It 
is  believed  that  this  difference  is  largely  due  to  the  inadequacy  of  the  four 
elements  used  here  to  model  the  increasingly  contorted  displacement  and  stress 
fields  which  occur  with  increasing  crack  speed.  An  increase  in  the  number  of 
triangular  elements  in  the  angular  direction  might  be  expected  to  remedy  this 
deficiency.  The  inadequacy  of  the  element  refinement  is  further  evidenced 


bv  the  noticeable-  disturbance  in  the  values  accompanying  each  remeshing. 


v=0.6C  The  results  for  this  largest  considered  crack  speed  are  presented 
in  Figure  5.  Here  again  the  dashed  line  is  the  steady-state  infinite  strip 
solution  and  the  solid  curve  the  propagation-eigen-function  solution.  The 
solid  square  points  are  Che  values  from  the  stationary-eigen-function  ele- 

s 

ment  solution,  and  are  obtained  from  the  (solid  circular  points)  using 
equation  19.  The  solid  triangular  points  are  also  K^.  values  from  the  sta¬ 
tionary-eigen-function  solution  but  were  obtained  through  fitting  the  near 
field  nodal  displacements  ( )  with  the  propagation-eigen-functions  using 
the  method  described  previously.  Nineteen  eigen-functions  plus  the  one  rigid 
body  mode  were  used  in  this  fitting  procedure.  Both  of  the  above  procedures 
for  treating  the  stationary-eigen-function-solution  yield  results  which  agree 
quite  well  with  the  propagacion-eigen-function-solution  but  clearly  the  one 
represented  by  equation  19  is  preferred  due  to  the  ease  of  application. 

The  results  from  the  quarter-point  isoparametric  elements  are  plotted 
in  Figure  5  using  open  circular  points.  Unlike  the  results  presented  at 
v=0.4Cs  and  v-0.2Cs>  the  values  of  plotted  here  were  obtained  through  the 
J-integral  as  defined  by  (16).  Though  the  results  based  on  the  global  energy 
balance  procedure  are  not  presented  here,  they  were  found  to  agree  quite 
well  with  those  using  the  J-integral  except  that  there  was  substantially  more 
"noise".  This  "noise"  became  particularly  apparent  at  positions  where  remeshing 
was  required. 

From  these  calculations  involving  crack  propagation  speed  ranging  from 

0.  If.  to  O.ui:  ,  it  appears  that  for  these  crack-speed  ranges  the  propag.at  ion- 
s  s 

eigen—! utu:  t  ion  element  has  no  significant  advantages  over  the  stat ionarv-e L gen- 
:unct  ien  element  ei  tiler  in  terr  oi  iccurncv  or  in  terms  of  ease  of  computing 
kj..  other  than  being  more  theoretical  lv  consistent  and  appealing. 


Since  there 


is  additional  effort  in  dealing  with  the  non-symwe t r ic  matrices  which  accompany 
the  propagation-eigen-function  formulation,  it  seems  the  stationary-eigen-func¬ 
tion  element  might  be  more  effectively  implemented  to  obtain  results  within 
a  tolerable  engineering  accuracy.  The  large  number  of  eigen- functions  used  in 
these  special  elenent  formulations  in  order  to  produce  no  spurious  kinematic  de¬ 
formation  modes  is  believed  to  be  the  reason  for  the  high  degree  of  accuracy 
attained  with  the  stationary-eigen-f unction  element  in  the  analysis  of  fast 
crack  propagation.  It  was  observed  that  the  use  of  four  quarter-point  iso¬ 
parametric  elements  was  quite  adequate  at  v=0.2C  ,  but  that  by  v=0.4C  the 

s  s 

accuracy  had  become  marginal  (5-10%  difference  from  special  element  results). 

In  considering  the  effect  of  crack  speed  on  the  isoparametric  element  solution 
accuracy,  it  should  be  kept  in  mind  that  experimentally  measured  crack  speeds 
often  fall  below  0.30^.  If  one  needs  to  consider  higher  speeds,  mesh  refinement 
seems  to  be  necessary. 


Application  to  Prediction  of  Crack  Growth 
In  the  computations  of  [1,3,5]  and  in  those  of  the  previous  section,  the 
loading  of  the  body  and  the  crack  growth  history  are  used  as  input  to  the 
analysis.  The  output  of  each  of  these  "generation  phase"  computations  is 
as  a  function  of  time,  crack  length  or  crack  speed.  The  subject  of  this  sec¬ 
tion  are  computations  of  a  reverse  nature.  That  is,  the  input  to  the  computation 
is  the  loading  of  the  body  and  a  crack  growth  history.  This  type  of  analysis 
is  referred  to  as  an  "application  phase"  computation. 

One  feature  of  the  special  'igen-f unct ion  elements  which  has  proved  to 


be  quite  useful  in  the  application  type  analvsis  is  thar  in  addition  to 

3KI  32Kl 

being  computi*d  direct]'/  from  <i  one  also  has  — —  and  - tt  as  a  result  of  the 

Jt  3t2 

min imi cat  ion  of  the  error  functionals  of  equation  (10). 
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Having  these  quantities  at  some  time,  say  t=t^,  simplifies  the  prediction  of 
crack  velocity  for  the  subsequent  time  step  since  the  value  of  can  extrap¬ 
olated  using  a  Taylor  series  expansion.  More  precisely,  in  order  to  determine 
AI=Z2_“^  using  a  crack  growth  criterion  based  on  one  needs  to  have  an 
average  K^.  for  the  interval  (t^.t^).  For  the  propagation-eigen- function 
element,  this  average  is  predicted  in  terms  of  quantities  at  t^  by: 

KiP  ■  +  t  8?<V +  ^  <21> 

3Kr  •  ••  32fCi 

where  3? = -  and  B?  =  - Having  K,.  one  uses  the  criterion  to  obtain  v 

1  3t  1  3  c  2  Ip 

which  is  then  used  to  obtain  AZ  =  vAt.  Having  AZ  the  finite  element  mesh 
for  the  crack  length  can  be  generated  and  the  solution  at  t^  obtained.  It 
should  be  noted  that  a  similar  prediction  procedure  to  (21)  exists  for  the 
stationary  eigen-function  element.  Through  the  differentiation  of  (19)  with 
respect  to  time  one  can  obtain: 

KiP '  6!ai> +  f  +  ^  <22> 

An  application  phase  calculation  for  a  double  cantilever  beam  specimen, 

identical  to  specimen  No.  U  of  Kalthoft',  et  al  [6],  has  been  completed  using 

the  K  versus  £  relation  show  in  Figure  6  as  a  crack  growth  criterion.  This 

curve  and  the  crack  initiation  fracture  toughness  for  the  computation  (K^  = 
-3/2 

2.32  MNm  “)  arc  identical  to  the  ones  used  by  Kobayashi  [7].  The  predicted 

kj.  $s  t  ,  Z  vs  t  and  Z  (or  v)  vs  t  which  were  obtained  using  the  propagation- 

eigen-funct ion  element  are  shown  in  Figure  7  along  with  the  corresponding 

experimental  results  of  Kalthoff,  et  aL  [6]  and  numerical  results  of  Kobayashi 

2  3 

(71.  lie  present  plane  stress  analysis  used  l".=3330MN/m  ,  v=0.33  and  .:  =  1172kg/m  . 


it  can  he  seen  in  Figure  7  that  the  ■  vs  t  results  are  virtuallv  indistin- 


_  n 


guishable  from  the  experimental  results  and  that  even  the  v  vs  t  results  agree 
quite  well.  In  comparing  the  curves  it  is  seen  that  there  is  quite  good 
agreement  for  approximately  the  first  half  of  the  total  crack  growth.  At 
about  midway,  however,  the  computed  increases  briefly  while  the  experimental 
values  drop.  Before  arrest  occurs,  there  is  again  good  agreement,  however. 

One  possible  explanation  for  the  above  disparity  is  that  the  Araldite 

B  used  by  Kalthoff  shows  some  rate  dependence  even  though  it  was  selected 

over  similar  materials  because  of  its  relatively  low  rate  dependence  [6]. 

For  example,  the  dynamic  elastic  constants  quoted  in  [6]  are  E=3660  and  v=0.39, 

2 

whereas  the  static  values,  used  in  the  analysis,  are  E=3380  MN/m  and  v=0.33. 
While  this  rate  dependence  must  surely  have  some  effect  on  the  specimens  res¬ 
ponse,  it  has  been  observed  through  numerical  experiments  that  other  aspects 
of  the  experimental  procedure  and  numerical  modeling  procedure  also  have  quite 
large  effects.  The  results  of  these  sensitivity  studies  are  being  presented  in 
a  companion  paper  [13]. 

CONCUR  TON 

The  comparison  of  the  propagation-eigen-function  element,  the  stationary- 
eigen-function  element,  and  the  quarter-point  isoparametric  element  in  terms  of 
accuracy  and  efficiency  leads  to  several  conclusions.  The  two  eigen-func¬ 
tion  elements  showed  similar  accuracy  at  all  crack  speeds  up  to  0.60^  while 
the  isoparametric  element  model  started  showing  significant  differences  be¬ 
tween  0.2C  and  0.4C  .  While  the  quarter-point  isoparametric  elements  were 
the  Least  expensive  to  use  (  2/1  the  cost  of  the  sta t ionarv-e igon-f unc t ion 
element  and  1/2  that  of  the  prop.ig.it  i  on-e  i  gen- f  une  t  ion  element)  their  sen- 
siti.’itv  to  crack  speed  and  the  need  to  use  indirect  methods  to  obtain  N 
reduce  their  initi.it  a  1 1  :act  i  veno  as  .  Furthermore,  for  application  tvpe 
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analyses  the  Lack  of  and  k^  may  complicate  the  use  of  isoparametric 
type  elements  to  the  point  that  thev  lose  their  cost  advantage.  Another  con¬ 
clusion  to  be  drawn  from  the  comparison  of  the  models  is  that  the  stationary- 
eigen-function  element  has  all  the  attractive  features  of  the  propagation 
eigen-function  element  without  the  disadvantage  of  non-symmetr ic  stiffness 
matrices.  Also,  because  the  matrices  are  not  crack  speed  dependent,  they 
do  not  need  to  be  recomputed  each  time  crack  speed  changes.  Finally,  the 

eigen-function  elements  seem  particularly  attractive  for  application  phase 

3Kj  3 2K 

analysis  since  at  each  time  step  one  has  not  only  K  but  also  ~r—  and  — — r-; 

1  dt  o  ^  Z. 

thus  simplifying  the  prediction  of  crack  behavior  in  the  subsequent  time  step. 

The  application  phase  anlaysis  of  the  double  cantilever  beam  specimen 
illustrated  the  utility  of  the  propagation-eigen-function  element  for  pre¬ 
dicting  crack  growth  behavior  using  versus  v  as  a  crack  growth  criterion. 
Since  the  stationarv-eigen-function  element  showed  very  good  agreement  with 
the  propagat ion-eigen- function  element  in  the  generation  type  analyses  and 
since  the  predictive  logic  for  application  type  analyses  is  the  same  for  both 
element  types,  it  seems  the  stationary-eigen-function  element  is  also  well 
suited  to  application  computations. 

The  propagation-eigen-function  element,  even  if  more  expensive  to  use  than 
the  other  two  special  elements  discussed  above,  is  nevertheless  more  consistent 
and  appealing  in  terms  of  theoretical  formulation,  and  application  to  basic 
research  in  dynamic  crack  propagation  in  finite  bodies.  Due  to  this  reason,  ex¬ 
tensive  dynamic  fracture  studies,  of  both  'generation'  and  'application'  type, 
on  laboratory  specimens  such  as  the  rectangular  double  cantilever  beam  (RDCD) 
tapered  double  cantilever  beam  (TDCB) ,  and  edge  crack  specimen,  were  conducted 
bv  the  authors  using  the  presently  desoribed  propngation-eigon-funecion  special 


.V. 


element . 


These  numerical  results  were  compared  with  available  experimental  data. 

A  detailed  presentation  of  these  results  is  made  in  a  companion  paper  [13],  in 
which  the  effects  of  specimen  geometry,  input  crack-velocity  history,  and  in¬ 
put  dynamic  fracture  toughness  property  data,  are  discussed  in  detail. 


acknowledgements 


This  work  was  supported  by  the  U.S.O.N.R.  under  contract  No.  N00014-78-7636 . 
The  authors  are  grateful  to  Dr.  N.  Perrone  for  his  timely  encouragement.  The 
authors  also  thank  Ms.  M.  Eiteman  for  her  care  in  the  preparation  of  this 
manuscript . 

REFERENCES 

[1]  Atluri,  S.N.,  Nishioka,  T.,  and  N'akagaki,  M. ,  "Numerical  Modeling  of 
Dynamic  and  Nonlinear  Crack  Propagation  in  Finite  Bodies,  by  Moving 
Singular  Elements"  In  Nonlinear  and  Dynamic  Fracture  Mechanics,  Ed.  N. 
Perrone,  and  S.N.  Atluri  AMD-Vol  35,  ASME,  NY,  1979,  pp  37-66. 

[2]  Nishioka,  T.,  and  Atluri,  S.N.,  "Numerical  Modeling  of  Dynamic  Crack 
Propagation  in  Finite  Bodies,  by  Moving  Singular  Elements,  Part  I  - 
Theory",  Journal  of  Applied  Mechanics,  ASME,  1980  (In  Press). 

[3]  Nishioka,  T.,  and  Atluri,  S.N.,  "Numerical  Modeling  of  Dynamic  Crack 
Propagation  in  Finite  Bodies,  by  Moving  Singular  Elements,  Part  II  - 
Results"  Journal  of  Applied  Mechanics,  ASME,  1980  (In  Press). 

[4]  Nilsson,  F.,  "Dynamic  Stress- Lntens i ty  Factors  for  Finite  Strip  Prob¬ 
lems",  Int.  Jnl.  of  Fracture,  Vol.,  No.  4,  1972,  pp  403-411. 

[5]  Nishioka,  T.,  and  Atluri,  S.N.,  "Efficient  Computational  Techniques 
for  the  Analysis  of  Some  Problems  of  Fracture  in  Pressure  Vessel 
and  Piping"  to  be  presented  at  ASME  Pressure  Vessels  and  Piping 
Conference,  San  Francisco,  Aug.  1980. 

[6|  Kalthoff,  J.F.,  Be  inert,  J.,  Winkler,  S.,  "Measurements  of  Dynamic  Stress 
Intensity  Factors  for  Fast  Running  and  Arresting  Cracks  in  Double-Canti- 
lever-Be.im  Specimens",  Fast  Fracture  and  Crack  Arrest,  ASTM  STP  627,  C.T. 
Haim  and  M.F.  Kanninen,  Eds.  ASTM,  1977,  pp  IM-176. 

[  7  J  Soil  sea  sit  i ,  A.S.,  "Dynamic  Fracture  Analysis  By  Dynamic  Finite  Element 

Method-Oenerat ion  and  Propagation  Analvses",  In  Nonlinear  and _ Dynamic 

Fracture  Meehan i cs,  Ed.  N.  Perrone.  and  S.N.  Atluri,  AMD-Vol  35,  ASME, 

NY,  lv7ri,  pp  IJ-j,,. 


[S]  Hcnshell,  R.D.,  Shaw,  K,C.,  "Crack  Tip  Finite  Elements  are  Unnecessary", 

Int.  J.  Num.  Meth.  Engnq.,  Vol  9,  495-507,  1975. 

[9]  Barsoum,  R.S.,  "On  the  Use  of  Isoparametric  Finite  Elements  in  Linear 
Fracture  Mechanics",  Int.  .J .  Num.  Moth.  Rngng.,  Vol  10,  25-37,  1970. 

[10]  Kishimoto,  K. ,  Aoki,  S.,  Sakata,  M.,  "Dynamic  Stress  Intensity  Factors 
Using  J-Integral  and  Finite  Element  Method”,  Engng.  Fracture  Mech.,  Vol 
13,  pp  387-394  (1980). 

[11]  Pian,  T.H.H.,  Tong,  P.,  "Basis  of  Finite  Element  Methods  for  Solid  Continua", 
Int.  J.  Num.  Heth.  Engng.,  Vol  1,  3-28  (1969). 

[1-1  Aciienbach,  J.D.  and  Bazant,  Z.P.,  "Elastodynaraic  Near-Tip  Stress  and 
Displacement  Fields  for  Rapidly  Propagating  Cracks  in  Orthotropic 
Materials",  J.  Appl.  Mech.,  42,  p  183  (1975). 


[13] 


Nishioka, 

Generation 

Mechanics, 


T.,  Atluri ,  S.N.,  "Numerical  Analysis  of  Dynamic  Crack  Propagati 
&  Prediction  Studies"  submitted  to  Jnl.  of  Engineering  Fracture 


on: 


Figure  1.  Finite  element  mesh  for  the  numerical  approximation  to  the  infinite 
strip  problem  of  Nilsson  [4], 

Figure  2.  Illustration  of  the  shifting/remeshing  procedure  for  modeling  crack 
growth  . 

Figure  3.  Nondimensionalized  K  vs  7  for  constant  velocity  crack  propagation 
(v=0.2Cs)  in  the  model  of  Figure  1. 

Figure  4.  Nondimensionalized  K^.  vs  7  for  constant  velocity  crack  propagation 
(v=0.4Cs)  in  the  model  of  Figure  1. 

Figure  5.  Nondimensionalized  vs  E  for  constant  velocity  crack  propagation 
(v=0.6Cs)  in  the  model  of  Figure  1. 

Figure  6.  vs  t  used  as  the  crack  growth  criterion  for  the  application 

(propagation)  phase  calculation  of  Figure  7. 

Figure  7.  Application  phase  analysis  of  the  double  cantilever  beam  specimen 
(No.  4)  of  Kalthoff  et  al  [6]. 
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